#Replication Code for
#Hansen, Eric R., and Andrew Tyner. "Educational Attainment and Social Norms of Voting." Forthcoming at Political Behavior.
#Compiled by Eric Hansen, 9/10/19 (ehansen4@luc.edu)
#Using R version 3.4.3, running on Windows 7 x64


rm(list=ls(all=TRUE))
set.seed(4198347)
#install.packages("Matching")
#install.packages("GenMatching")
#install.packages("rgenoud")
#install.packages("mediation")
#install.packages("cobalt")
#sessionInfo()

library(ggplot2)
library(Matching)
library(rgenoud)
library(mediation)
library(foreign)
library(cobalt)

####code we used to match variables from ANES 2016 Time Series data set to the validated vote data. Time Series variables are the same as defined in the code used for Study 1 in the paper. We've cut out these steps for replicators by providing the merged data.

#data <- read.csv("C://Users/ehansen4/Dropbox/Eric and Tyner/Class and Motivation/ANES 2016/anes_2016_voteval_merge.csv")
#data <- data.frame(data$valvote, data$college, data$age, data$age2, data$female, data$white, data$black, data$asian, data$hispanic, data$foreignborn, data$religattend, data$pidstr, data$democrat, data$gop, data$duty, data$scispend)
#data <- na.omit(data)
#write.csv(data, file="C://Users/ehansen4/Desktop/matchdata.csv")

data<- read.csv("C://Users/ehansen4/Desktop/matchdata.csv")
data$X<-NULL

#matched covariates
attach(data)
X <- (cbind(data.age, data.age2, data.female, data.white, data.black, data.asian, data.hispanic, data.foreignborn, data.religattend, data.pidstr, data.democrat, data.gop))

BalanceMat <- (cbind(data.age, data.age2, data.female, data.white, data.black, data.asian, data.hispanic, data.foreignborn, data.religattend, data.pidstr, data.democrat, data.gop))

#fyi, this might take a while. runtime on my Windows 7 desktop PC was about half an hour.
genout <- GenMatch(Tr=data.college, X=X, BalanceMatrix=BalanceMat,M=1,
                   pop.size=1000, max.generations=50, wait.generations=1)

mgen1 <- Match(Y=data.valvote, Tr=data.college, X=X, Weight.matrix=genout, ties=F)
summary(mgen1)

#balance check
new.names <- data.frame(old = c("data.age", "data.age2", "data.female", "data.white", "data.black", "data.asian", "data.hispanic", "data.foreignborn", "data.religattend", "data.pidstr", "data.democrat", "data.gop"),
    new = c("Age", "Age X Age", "Female", "White", "Black", "Asian", "Hispanic", "Foreign Born", "Religious Attendance", "PID Strength", "Democrat", "Republican"))

love.plot(bal.tab(mgen1, data.college~data.age + data.age2 + data.female + data.white + data.black + data.asian + data.hispanic + data.foreignborn + data.religattend + data.pidstr + data.democrat + data.gop, data=data),
          title=" ", var.names=new.names, shapes = c(16, 22), colors = c("black", "black"))


detach(data)


data$rowID<-as.numeric(rownames(data))
matched.data <- data.frame(mgen1$mdata,mgen1$index.treated,mgen1$index.control)

#causal mediation model with duty
matched.data$TreatDuty <- data$data.duty[match(matched.data$mgen1.index.treated, data$rowID)]
matched.data$ControlDuty <- data$data.duty[match(matched.data$mgen1.index.control, data$rowID)]
matched.data$duty <- ifelse(matched.data$Tr == 1, matched.data$TreatDuty, matched.data$ControlDuty)

model.m <- lm(duty ~ Tr, data = matched.data)
model.y <- glm(Y ~ Tr + duty, data = matched.data, family = binomial(link="probit"))
summary(model.m)
summary(model.y)

acme <- mediate(model.m, model.y, sims=1000, treat="Tr", mediator="duty")

summary(acme)
plot(acme)

sens.out <- medsens(acme, rho.by = 0.1, effect.type = "both", sims = 100)
summary(sens.out)

plot(sens.out, sens.par = "rho", r.type = "total")
plot(sens.out, sens.par = "R2", r.type = "total", sign.prod="positive")

#placebo model with science spending

matched.data$TreatSci <- data$data.scispend[match(matched.data$mgen1.index.treated, data$rowID)]
matched.data$ControlSci <- data$data.scispend[match(matched.data$mgen1.index.control, data$rowID)]
matched.data$scispend <- ifelse(matched.data$Tr == 1, matched.data$TreatSci, matched.data$ControlSci)

model.pm <- lm(scispend ~ Tr, data = matched.data)
model.py <- glm(Y ~ Tr + scispend, data = matched.data, family = binomial(link="probit"))
summary(model.pm)
summary(model.py)


p.acme <- mediate(model.pm, model.py, sims=1000, treat="Tr", mediator="scispend")

summary(p.acme)
plot(p.acme)

p.sens.out <- medsens(p.acme, rho.by = 0.1, effect.type = "both", sims = 100)
summary(p.sens.out)


plot(p.sens.out, sens.par = "rho", r.type = "total")
plot(p.sens.out, sens.par = "R2", r.type = "total", sign.prod="positive")

